Transport anisotropy of the pnictides studied via 
Monte Carlo simulations of a Spin-Fermion model 
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An undoped three-orbital Spin-Fermion model for the Fe-based superconductors is studied via 
Monte Carlo techniques in two-dimensional clusters. At low temperatures, the magnetic and one- 
particle spectral properties are in agreement with neutron and photoemission experiments. Our main 
results are the resistance vs. temperature curves that display the same features observed in BaFe2 As2 
detwinned single crystals (under uniaxial stress), including a low-temperature anisotropy between 
the two directions followed by a peak at the magnetic ordering temperature, that qualitatively 
appears related to short-range spin order and concomitant Fermi Surface orbital order. 

PACS numbers: 



Introduction. In early studies of Fe-based supercon- 
ductors [1], it was widely assumed that Fermi Surface 
(FS) nesting was sufficient to understand the undoped- 
compounds magnetic order with wavevector Q — (7r,0) [2] 
and the pairing tendencies upon doping. Neutron 
scattering reports of spin-incommensurate order |3] are 
in fact compatible with the nesting scenario. How- 
ever, several recent experimental results cannot be ex- 
plained by FS nesting, including (i) electronic "ne- 
matic" tendencies in Ca(Fei_a;COa;)2As2 [1]; (ii) orbital- 
independent superconducting gaps [5] found using laser 
angle-resolved photoemission (ARPES) spectroscopy on 
BaFe2(Aso.65Po.35)2 and Bao.6Ko.4Fe2As2; and, more im- 
portantly for the investigations reported here, (in) the 
report of local moments at room temperature (T) via 
Fe X-ray emission spectroscopy [5J [7] . Considering these 
experiments and others, a better characterization of the 
pnictides is that they are in the "middle" between the 
weak and strong Coulomb correlation limits [SUlOj. Be- 
cause this intermediate Hubbard U range is difficult for 
analytical approaches, there is interest in the develop- 
ment of simpler models that can be studied via compu- 
tational techniques to provide insight into such a difficult 
coupling regime. The Hartree-Fock (HE) approximation 
to the Hubbard model [H] cannot be applied at room 
T since HF approximations only lead to non-interacting 
fermions above the ordering temperature Tn, and thus 
the local moment physics [S] cannot be reproduced [T^]. 

Recently, a Spin-Fermion (SF) model for the pnictides 
has been independently proposed by Lv et al. [I3] and 
Yin et al. |14j . The model, very similar to those widely 
discussed for manganites, originally involved itinerant 
electrons in the xz and yz d-orbitals coupled, via an on- 
site Hund interaction, to local spins (assumed classical) 
that represent the magnetic moment of the rest of the Fe 
orbitals (considered localized). The Hund interaction is 
supplemented by a nearest-neighbor (NN) and next-NN 



(NNN) classical Heisenberg spin-spin interaction. This 
SF model has interesting features that makes it quali- 
tatively suitable for the pnictides, particularly since by 
construction the model has itinerant electrons in interac- 
tion with local moments [51 [7] at all temperatures. 

Phenomenological SF models have been proposed be- 
fore for underdoped cuprates, with itinerant fermions 
representing carriers locally coupled to classical spins rep- 
resenting the antiferromagnetic order parameter. These 
investigations unveiled stripe tendencies [T^, ARPES 
and optical conductivity results [T6| similar to experi- 
ments, and even the dominance of the d-wave channel in 
pairing |17) . Thus, it is natural to apply now these ideas 
to the Fe superconductors. 

As remarked already, SF models are also mathemat- 
ically similar to models used for the manganites |18j . 
Then, all the experience accumulated in the study of Mn- 
oxides can be transferred to the analysis of SF models 
for Fe-superconductors. In particular, one of our main 
objectives is to study for the first time a SF model for 
pnictides employing Monte Carlo (MC) techniques, al- 
lowing for an unbiased analysis of its properties. More- 
over, to test the model, challenging experimental re- 
sults will be addressed. It is known that for detwinned 
Ba(Fei_a;Co2:)2As2 single crystals, a puzzling transport 
anisotropy has been discovered between the ferromag- 
netic (EM) and antiferromagnetic (AFM) directions |19) . 
In addition, the resistivity vs. T curves display an un- 
expected peak at Tn~130 K, and the presumably weak 
effect of an applied uniaxial stress [19] still causes the 
anisotropy to persist well beyond T^. However, recent 
neutron results suggest that the transport anisotropy 
may be actually caused by strain effects that induce a 
shift upwards of the tetragonal-orthorhombic and Tn 
transitions [20l [21], as opposed to a spontaneous ro- 
tational symmetry-breaking state not induced by mag- 
netism or lattice effects. Then, theoretical guidance is 
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needed. While the low-T anisotropy was already ex- 
plained as caused by the coupling between the spins and 
orbitals in the Q=(7r, 0) state |22], the full transport 
curves at finite T define a challenge that will be here 
addressed for the nontrivial undoped case. 

Model and Method. The SF model [H [H] is given by 



pp 



Hund 



Hcis • 



(1) 



The first term -ffnopp describes the Fe-Fe hopping of itin- 
erant electrons. To better reproduce the band struc- 
ture of pnictides [Ij, three d-orbitals (ccz, yz, xy) will 
be used instead of two. The full expression for -ffnopp 
is cumbersome to reproduce here, but it is sufficient for 
the readers to consult Eqs.(l-3) and Table I of Ref. 
for the mathematical form and the actual values of 
the hoppings in eV's. The density of relevance used 
here is n=4/3 [531. The Hund interaction is simply 
-f^Hund=— >/h X^ia ' '^'^h Si the classical spin at site 
i (|Si|=l), and s"i the itinerant-fermion spin of orbital 
a |24| . The last term ifneis contains the spin-spin interac- 
tion among the localized spins i?Hcis=>/NN ' + 
•^NNN E((im» Si • where () ((())) denotes NN (NNN) 
couplings. The particular ratio Jnnn / '/nn=2/3 was used 
in all the results below, leading to (tt, 0)/(0, tt) mag- 
netism [25]. Any other ratio Jnnn/^/nn larger than 1/2 
would have been equally suitable for our purposes. 

The well-known Monte Carlo (MC) technique for SF 
models [iH] will be here used to study i?sF at any T. 
In this technique, the acceptance-rejection MC steps are 
carried out in the classical spins, while at each step a 
full diagonalization of the fermionic sector (hopping plus 
on-site Hund terms) for fixed classical spins is performed 
via library subroutines in order to calculate the energy 
of that spin configuration. These frequent diagonaliza- 
tions render the technique CPU-time demanding. The 
simulation is run on a finite 8x8 cluster with periodic 
boundary conditions (PBC) and uses the full i?sF model 
for the MC time evolution and generation of equilibrated 
configurations for the classical spins at a fixed T [55]. 
However, for the MC measurements those equilibrated 
configurations are assumed replicated in space but dif- 
fering by a phase factor such that a better resolution 
is reached with regards to the momentum k. Since a 
larger lattice with more eigenstates gives a more con- 
tinuous distribution of eigenenergies, the procedure then 
reduces finite-size effects in the measurements. This well- 
known method is often referred to as "Twisted" Bound- 
ary Conditions (TEC) 27 . In practice, phases $ are 
added to the hopping amplitudes, schematically denoted 
as "t", at the boundary via iTSC=e**t, with ^—2-Km/M 
(to=0,1,...,A/ — 1) such that the number of possible mo- 
menta in the x or y directions becomes L=8xM. 

Results. Figure [ija) contains MC results for the 
structure factor of the classical spins, defined as 
5(q) = ^Ei j(Si • Sj)e^i-('-j) (N = number of sites), il- 
lustrating the development of Q=(7r, 0) magnetic order 
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FIG. 1: (color online) (a) Classical spins structure factor 
S(7r,0) (normalized to 1) vs. T, for the Jh's indicated, us- 
ing the PBC 8x8 cluster and Jnn=0.015. The oscillations 
in the data are indicative of the error bars, (b) Density of 
states N{(jj) of each orbital (/i= chemical potential), using 
TBC with L=512, at T^O K and Jh=0.1 eV, for the perfect 
(tTjO) magnetic state, (c) Resistance R vs. L (TBC 8x8) for 
the FM and AFM directions of the perfect (tt, 0) magnetic 
state (Jh=0.1 eV). (d) The occupation at the FS n(/i) (see 
text) of the three orbitals vs. T, using L=256. A coupling 
Jnn^O.OIG (0.014) along the x (y) axis was used (see text). 

as T is reduced. Since a ratio Jnnn />/nn > 1 /2 is "frus- 
trating" , finding Q-order at Jh=0 is not surprising, but 
Fig. [T] shows that this order remains stable turning on Jjj 
in the range investigated, as opposed to inducing transi- 
tions to other states. The chosen value of Jnn in Fig.[l](a) 
leads to a Tn similar to that in BaFe2As2. The low- 
T orbitally-resolved electronic density-of-states (DOS) is 
in Fig. [ijb). The Q magnetic order opens a pseudogap 
(PC) in the yz orbital, while the others are not much 
affected. This PC generation was previously discussed 
when contrasting theory and ARPES experiments [21] 
and it should not be confused with long-range orbital- 
order, that in this SF model occurs at Jh~0.4 or larger. 

Figure [ijc) contains the evolution of the 8x8-cluster 
resistance R increasing the number of momenta via the 
TBC, calculated via standard procedures [T5J[55]. While 
the ratio of i?'s along the FM and AFM directions is al- 
ways > 1, i.e. qualitatively correct, TBC with L=256 
is needed to reach stable values. In addition, the occu- 
pation of the three orbitals at the FS (Fig. W[d)) was 
defined as n{n)= J den(e)/3e''('-'') /(I + e'3(^-M)p^ involv- 
ing the ^-derivative of the fermionic population. As T 
increases and the Q order weakens, the xz-yz orbitals 
populations converge to the same values. 

The results of Fig. [T] and others below, were obtained 
introducing a "small" explicit asymmetry along the x and 
y axes for Jnn, namely a generalized {a=x,y) was 
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FIG. 2: (color online) A{k,ui) at (TBC 8x8 L=512). 

The model used (see Ref. [23l includes the staggered As mod- 
ulation out of the FeAs layer. Thus, our results are in the 
folded Brillouin zone convention, and for this reason two elec- 
tron pockets (as opposed to just one) are centered at X in 
the panels above. The pocket elongated vertically at X would 
corresponds to a pocket at Y in the unfolded convention if 
the As modulation is considered via a quasicrystal momen- 
tum [23]. (a-c) are for Jh=0.1 eV. Red, green, and blue are 
for the xz, yz, and xy orbitals, respectively, (a) T='10K, be- 
low Tn. The (tt, 0) magnetic order induces j/z-orbital electron 
satellite pockets, (b) T^llOKr^Tj^. In this regime, the MC 
configurations display small coexisting patches of (tt, 0) and 
(0,7r) order (see text), creating almost symmetric xz and yz 
features around P. (c) Large T=360 K, with no remnants of 
the (tt, 0) order. The FS becomes a broaden version of the 
non-interacting FS at Jh=0, shown in (d) at T^lOO K. 

used. Its purpose is to mimic the orthorhombic distortion 
and strain effects [2 [^Ht 2d^ and judge if the present cal- 
culations reproduce transport experiments [I9H2T] . Con- 
sider the ratio '''nn='^nn/'^nn- Using the dependence of 
the hopping amplitudes with the distance u between el- 
and p-orbitals i.e. tpd^l/u^^^ 13D], the angles involved 
in the Fe-As-Fe bonds, the low-T lattice parameters 0, 
fourth-order perturbation in the hoppings for the Fe-Fe 
superexchange, and, more importantly, neglecting con- 
tributions of the yz orbital that is suppressed at the 
Fermi level |2S] as long as (tt, 0) spin fluctuations domi- 
nate leads to a crude estimation Since this is 
likely an upper bound, the ratio used in our MC simula- 
tion rNN^l-14, assumed to be temperature independent, 
is reasonable. Other crude estimations including the di- 
rect Fe-Fe hoppings |30) or employing the lattice param- 
eters under pressure |20) lead to similar ratios. Then, our 
asymmetry value is qualitatively realistic. 

Previous investigations [TT] showed that the T=0 HF 
approximation to the undoped multiorbital Hubbard 
model can reproduce neutron diffraction results and 
ARPES data. A similar degree of accuracy should be 
expected from any reasonable model for the pnictides, 
including -ffsF- To test this assumption, the one-particle 
spectral function yl(k,w) was calculated, and the FS at 
different T's is shown in Fig. [2j contrasted against the 
low-T fermionic non-interacting limit Jn—O. At low-T 
in the ordered state the expected asymmetry between 
the (tt, 0) and (0, tt) electron pockets is observed (not 
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FIG. 3: (color online) Resistance R of the SF model calculated 
via MC simulations, at Jh=0.1 eV and using the L—256 TBC 
8x8 cluster. Tn is indicated, and the magnetic susceptibili- 
ties xs are also shown. The FM-AFM directions asymmetry 
is evident at low T (note that the FM and AFM labels sim- 
ply refer to the y and x directions, respectively, and not to 
fully FM or AFM spin configurations). As T increases the 
symmetry is restored, and the curves display a peak at Tn. 
A small symmetry-breaking difference rNN=1.14 is used (see 
text). In green (dashed line) are the results for random spin 
configurations, showing that their R is smaller than in the MC 
results near Tn. The width of the x peaks extend to ~1.5Tn, 
in agreement with neutron scattering for CaFe2As2 |32| . 

shown), and more importantly satellite pockets (with 
electron character) develop close to the F hole pockets 
(Fig.|2];a)), as in ARPES experiments ^,WP. Thus, the 
SF model studied here passes the low-T ARPES test. As 
T increases, at or well above Tn (Figs.[2jb) and (c)) the 
xz and yz differences are reduced and rotational invari- 
ance is recovered, albeit with a FS broader than in the 
non-interacting low-T limit (Fig. [2]jd)). 

R vs. T curves. Our most important result is the T 
dependence of R along the two axes, shown in Fig. [3] It 
is visually clear that these results are similar to the trans- 
port data of Ref. [TF', particularly after realizing that lat- 
tice effects, that cause the continuous raise of R with T in 
the experiments, are not incorporated in the SF model. 
A clear difference exists between the FM and AFM direc- 
tions at low T, induced by the (tt, 0) magnetic order that 
breaks spontaneously rotational invariance. At low T, 
this difference was understood in the Hubbard-model HF 
approximation [22j based on the reduction of the yz or- 
bital population (Fig. [l|d)). This explanation is equally 
valid in the SF model, and at low T the SF model and the 
Hubbard model, when treated via the HF approximation, 
lead to similar physics. 

The most interesting result in Fig. [3] is the develop- 
ment of a peak at Tn, and the subsequent slow conver- 
gence of R toward similar values along both directions 
with further increasing T (as already discussed, to model 
better the effect of uniaxial stress fTSJ , a weak symmetry- 
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FIG. 4: (color online) Resistance R calculated using a L=64 
TBC 16x16 cluster, with Jnn=0.032 (0.028) along the x (y) 
axis. The spin configurations are generated via Hucis, while 
the R measurements use the full Hsf, at Jh=0.1 eV. The rest 
of the notation is as in Fig. [3] 

breaking difference between the NN Heisenberg couplings 
along X and y was included). To our knowledge, this is 
the first time that the full R vs. T curve is successfully 
reproduced via computational studies. 

A study as in Fig. |3] using a larger cluster, e.g. 16 x 16, 
is not practical since the computer time grows like TV^ 
(A^=number of sites), leading to an effort 256 times 
larger. However, results as in Fig. [ija) indicate that 
the classical spins configurations generated merely by 
the spin-spin interaction could be qualitatively similar 
to those generated by the full SF model, as long as Jh 
does not push the system into a new phase. Thus, the 
MC evolution could be carried out with i?Hois only, while 
measurements can still be performed using the full diag- 
onalization of i?sF- Such measurements (very CPU-time 
consuming) must be sufficiently sparsed in the MC evolu- 
tion to render the process practical. This procedure was 
implemented on a TBC 16x16 cluster, with L=64 [33] . 
The results for R are in Fig. [4j and they show a remark- 
able similarity with Fig. [3j and with experiments. Thus, 
the essence of the R vs. T curves is captured by electrons 
moving in the spin configurations generated by -ffnois- 
Size effects are small in the range analyzed here. 

What causes the increase of R upon cooling before Tn 
is reached, displaying insulating characteristics? Since 
our results are similar to experiments, an analysis of the 
MC-equilibrated configurations may provide qualitative 
insight into their origin. In Fig. [sja), a typical MC con- 
figuration of classical spins is shown. The colors at the 
links illustrate the relative orientation of the two spins 
at the ends. The (tt, 0) long-range order is lost, but in- 
dividual spins are not randomly oriented. In fact, the 
state contains small regions resembling locally either a 
(7r,0) or (0,7r) order (short-range spin order), and S'(q) 
still displays broad peaks at those two wavevectors. In 
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FIG. 5: (color online) (a) Typical MC-generated classical spin 
configuration on an 8x8 cluster at Tn^IIO K, Jh=0.1 eV, and 
Jnn=0.016 (0.014) along the x (y) axis. The red (blue) lines 
denote AFM (FM) NN correlations, of intensity proportional 
to the width, (b) The dominant orbital at the FS at each site 
for the configuration used in (a), calculated using n{fi) as in 
Fig. 1(d). Red (green) denote the xz (yz) orbital, with a size 
proportional to the density. 

standard mean-field approximations there are no precur- 
sors of the magnetic order above Tn, but in the SF model 
there are short-range fiuctuations in the same regime. 

States with "spin patches" as in Fig. [5|a) lead to a 
concomitant patchy orbital order at the FS shown in 
Fig. [sjjb), where the most populated orbital at /i at 
every site, either xz or yz, are indicated. The orbital 
orientation suggests that the xz (yz) FS population fa- 
vors transport along the x (y) axis. The patchy order 
should have a resistance larger than that of a randomly- 
oriented spin background. This is confirmed by calcu- 
lating R vs. T in the absence of a guiding Hamiltonian, 
i.e. by generating random spin configurations. The re- 
sults are in Figs.[3]and|4]( green dashes) and their values 
are indeed below those of the peak resistance at Tn of 
the full SF model, i.e. with configurations as in Fig. [5] 
Then, the effect of strain coupled to short-range spin 
and FS orbital order appears to be the cause of the 
peak in the R vs. T curves [HH [35]. Using a smaller 
(but nonzero) anisotropy, the R-T curves display a con- 
comitantly smaller anisotropy, but still they have a small 
peak at Tn (not shown). Thus, the patchy states may 
also explain the insulating properties of Fei.osTe [M] and 
(Tl,K)Fe2-^Se2 [37] above Tn. 

Conclusions. The Spin-Fermion model for pnictides 
was here studied with MC techniques. The magnetic 
and ARPES properties of the undoped compounds are 
well reproduced. Including a small explicit symmetry 
breaking to account for strain effects, the resistance R 
vs. T curves are qualitatively similar to those observed 
for BaFe2As2 [IH], including a peak at Tn that appears 
caused by short-range spin and FS orbital order. In our 
calculations, the anisotropy above Tn exists only as long 
as a strain distortion exists, compatible with results for 
annealed BaFe2As2 samples [32 . This successful appli- 
cation of a SF model paves the way to more demanding 
efforts involving doped systems where anisotropy effects 
are stronger than in the undoped limit [19j . 
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